function [mode, c, cjac, user] =    confun_mean(mode, ncnln, n,nrowj, needc, w, cjac, nstate, user)
    
    %%%%%%%%%%%%%%%%%%%%%%%
    % INITIALIZE
    %%%%%%%%%%%%%%%%%%%%%%%

[~,~,~,~,~,n1,n2,k1,k2,~,~,p_grid_constraint,y_grid_constraint]=user_unpack(user);

sizew=(n1+k1)*(n2+k2);

    %%%%%%%%%%%%%%%%%%%%%%%
    % PREDICT SHAT 
    %%%%%%%%%%%%%%%%%%%%%%%

[shat_grid,S_grid,S_deriv1_grid,S_deriv2_grid]=predict_demand_share(p_grid_constraint,y_grid_constraint,w,user);

    %%%%%%%%%%%%%%%%%%%%%%%
    % COMPUTE CONSTRAINT
    %%%%%%%%%%%%%%%%%%%%%%%

c = S_deriv1_grid*w + (S_grid*w).*(S_deriv2_grid*w) - shat_grid + shat_grid.^2;

cjac = S_deriv1_grid + S_deriv2_grid.*(S_grid*w*ones(1,sizew)) ...
         + S_grid .* ((S_deriv2_grid*w)*ones(1,sizew)) ...
         - (S_grid) + 2*(shat_grid*ones(1,sizew)).*S_grid;       

    %
